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Theoretical alternatives to the Lotka-Volterra model 


The Lotka-Volterra model provides a phenomenological representation of bacterial growth and interactions, but our conclusions are not tied to this specific 
model. It is important to understand whether our central theoretical predictions generalize: how broadly do we expect a similar qualitative map of dynamical 
phases, where extinctions start to occur before the onset of fluctuations (as we increase either species pool size or interaction strength)? This empirically- 
observed ordering of phases is not self-evident: it is straightforward to construct few-species models that display fluctuations without extinctions, e.g. 
predator-prey pairs. Therefore, if our qualitative phase ordering appears across a range of many-species models and experiments, we may be seeing a 
broad mechanism, one that is presumably collective rather than driven by particular species. 


One interesting question is whether these collective dynamics emerge from many independent pairwise species interactions (as in our random Lotka- 
Volterra model), or whether they are driven by one or a few system-wide factors, such as public goods impacting all species. Ratzke et al. (74) performed 
bacterial experiments and introduced a different model where all interactions are mediated by modifications of the environmental pH by the bacteria, whose 
growth is in turn modified by pH. This pH-based model, and a variant, reproduced some experimental results in (14, 33). We slightly amend the model in (74) 
to represent continuous-time dilution and dispersal from the species pool (we have also simulated discrete daily dilutions and dispersal, with no impact on 
our conclusions). 
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Equation (2) represents the growth of bacterial abundance 
Wi 
, which follows a logistic equation with parameter 
modulated by the pH value 
a 
: growth is maximal when this value is equal to species /’s pH optimum 
bal 
, (drawn uniformly over [4.5,9.5]), with a tolerance given by 
p. = q 
. Inaddition, parameter 


(continuous-time equivalent to our 1:30 daily dilution) encompasses losses due to dilution, whereas 
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represents dispersal from the species pool. Equation (3) represents the change of pH induced by the bacteria, and the return toward the neutral pH value of 
7 due to dilution The distribution of 
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, drawn uniformly over [-c, c], thus determines the impact of bacteria on pH, and indirectly, the strength of interactions between bacteria. 


We find that this model reproduces some of the main predictions of the Lotka-Volterra model: as seen in Fig. S10, the phase diagram with its three phases 
is qualitatively preserved in this second model, which suggests that it is highly robust to variations in modelling assumptions. Nevertheless, we believe that 
the LV model better reproduces our experimental results, whereas the pH-based model was more adequate in (14, 33), plausibly due to different taxa 
(species that exhibit ecological suicide do not appear in our species pool here) and experimental conditions (lower dilution rate) that lead 
to stronger impacts of pH. Indeed, there are three points that make the LV model more plausible here: 


1) In our experiments, pH fluctuations are only moderately correlated (Pearson correlation coefficient: 0.54) with fluctuations in optical density or species 
composition (whereas the latter two are highly correlated), suggesting that pH is not the sole driving factor of stability here. (Fig S11A) 


2) The pH-based models in (14, 33) notably aimed to allow for ecological suicide, i.e. species growing then going extinctin monoculture, which does not 
occur in our experiment (In our initial species pool, we only retained species that survived in all of our growth conditions for the convenience of daily dispersal, 
Fig. S26). In the high-interaction strength regime and large initial species pool size regime, the pH-based model has two main outcomes: either many- 
species fluctuations, or total or near-total extinction with zero or few species surviving, and total biomass < 5% of carrying capacity (Fig S11B). Only the 
former behavior is seen in our experiments; the latter prediction of pH-based model was never observed in our experiments (Fig. 2C and Fig. S26). 


3) The LV model displays positive correlation between the number of surviving species and the intensity of fluctuations, as in our experiments (Fig. 4 and Fig. 
S11), whereas the pH-based model displays a negative correlation (we excludecommunities falling under the above extinction criterion, which were never 
observed in our experiment, Fig. S11C). The positive relation between diversity and fluctuations, validated in our experiments, is key to our theoretical 
argument. 


These three reasons lead us to retain the Lotka-Volterra model as our prime example (see Fig. S12 for a best-fit illustration). This cannot rule out a different 
pH-based model, or a combination of pH-based and direct interactions, as a good description of the biological mechanisms at play in our experimental 
setting. 

We reach two main conclusions: 

i) The random Lotka-Volterra model provides a straightforward prototype of a mechanism that fits our theoretical argument. The dynamical phases are 
unambiguously emergent, in the sense of being driven by species diversity and local pairwise interactions. The picture is as follows: the fluctuating 
phase must appear after extinctions and turnover in species composition, since fluctuations are driven by a tendency to jump between alternate sets of 
surviving species, each of themunstable. 

ii) We expect this central prediction to hold broadly under different modelling choices. Indeed, we find that the pH-based model displays an analogous set of 
regimes and transitions (Fig. S10), even though its predictions do not match some of our experimental results as accurately as the LV model (Fig. S11). 

In both models, increasing the number of species and the strength of their interactions, whether direct or mediated through an abiotic factor, will typically lead 
to extinctions before it leads to loss of stability. This results in a phase of partial coexistence preceding a phase of instability, no matter whether species 
diversity and extinctions are necessary (as in the random LV model) or not (as in the pH-based model) for the mechanism that drives fluctuations. 
We propose that this ordering, which is robust in many-species models but need not be in few-species models, can be an indicator of emergent collective 
behavior. 


Below is the code for our pH-mediated model: 


import numpy as np, pandas as pd, sys, matplotlib.figure as mpfig, matplotlib.pyplot as plt 
from scipy.integrate import solve_ivp 


EPS=10.**-8 ### Numerical threshold 
dispersal=10**-6 
def run( 
S=20, #Species number 
cp=100., #Coupling strength 
days=15, 
pHO=7, 
dilution=0.1, #if >0, daily dilution 
pc=2.5, #gate width 
kgrowth=10., #growth rate 
™kwargs 
): 
ppref=kwargs.get(‘ppref ,np.random.uniform(pHO-pc,pHO0+pc,S)) 
cs= kwargs.get('cs',np.random.uniform(-1,1,S)*cp ) 
x0=kwargs.get('x0',np.concatenate([[pHO],np.log(np.random.random(S)*0.01)}) ) 
def fun(p): 
ingate=(np.abs(p-ppref)<pc ).astype(‘float’) 
return ingate 
def eqs(t,x): 
p=x{0] 
x=np.exp(np.clip(x{1:],np.log(EPS),0)) 
dx= x * (1-x)* kgrowth*fun(p) 


if dilution<0: 
#Continuous-time dilution 
dxt= np.log(-dilution)*x 
dx+= dispersal 
dp=np.sum(cs*x) 
if dilution<0: 
dp+= np.log(-dilution)*(p-pHO) 
if p>14: 
dp=min(0,dp) 
if p<0: 
dp=max(0,dp) 
return np.concatenate([[dp],dx/x]) 


time=]] 

pH=[] 

Ns=]] 

for day in range(days): 
sol=solve_ivp(fun=eqs, t_span=(0,1.),yO=x0) 
pH=np.concatenate([pH,sol.y[0]]) 
if len(Ns): 

Ns=np.concatenate([Ns,np.exp(sol.y[1:].T)]) 
else: 

Ns=np.exp(sol.y[1:].T) 
time=np.concatenate([time,sol.ttday]) 
x0=sol.y[:,-1].copy() 
if dilution>0: 

#Discrete daily dilution 

x0[0]=(1-dilution)*pHO + dilution* x0[0] 

x0[1 :]=np.clip(xO[1:]+np.log(dilution),np.log(EPS),None) 

return time[:-1], H[:-1], Ns[:-1] 


def measure(ts,pH,Ns): 
Nf=Ns[-1] 
S=len(Nf) 
days=sorted(set(ts//1)) 


observation_period=int(max(days))//3 

Nmean=np.array([ np. mean(Ns[ts//1 ==z],axis=0) for z in days ]) 

nrel=Nmean/np.sum(Nmean,axis=1 ).reshape((-1,1)) 

pHmean=np.array([ np.mean(pH{ts//1 ==z]) for z in days ]) 

Ntot=np.sum(Nmean,axis=1) 

alive=np.where(Nmeanr][-1]>dispersal*S )[0] 

nalive=len(alive) 
CV=np.std(Nmear[-observation_period:,alive],axis=0)/(0.01+np.mean(Nmear[-observation_period:,alive],axis=0)) 
meanCV=np.mean(CV) 

fluct=(meanCV > 0.25) 


return {‘alive':float(nalive)/S, 'fluct’:float(fluct),} 


Ss=[1,2,3,4,6,8,10,12,15,18,21,24,30,36,42,48] 
cps=np.logspace(0,4,13) 
replicas=500 
days=20 
table=[] 
dilution=-0.033 
fname='pHdata.csv’ 
df=None 
for SinSs: 
for cp in cps: 
print 'RUN:, S, cp 
nrep=replicas 
if S==Ss[-1]: 
nrep=replicas*2 
for replica in range(nrep): 
print’ ‘replica 
opts={‘days':days,'cp':cp,'dilution':dilution,} 
ts, pH,Ns=run(S=S,**opts) 
dic={'S':S, 'cp':cp,'replica':replica} 
dic.update(measure(ts,pH,Ns)) 
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